README

In this R project, all the scripts for the manuscript Plankton blooms over the annual cycle shape trophic interactions under climate change are available.

All the analyses were done with R version 4.1.1 (2021-08-10) on the platform x86_64-apple-darwin17.0.

Note that Figure numbering might not be the same as the one in the manuscript.

R scripts

# 1. Download the data
source("Scripts/DownloadData.R")

# 2. Calculate the phenology metrics and the phyisical variables
source("Scripts/PhenologyZP.R")
source("Scripts/PhenologyPP.R")
source("Scripts/Physical.R")

# 3. Run the analysis
source("Scripts/PlanktonSuccession.R")
source("Scripts/Trends.R")
source("Scripts/GLM.R")

# 4. Run the supplementary analyses
source("Scripts/DifferenceStations.R")
source("Scripts/Nauplii.R")
source("Scripts/FreqPhytoProp.R")
source("Scripts/Validation_Trends.R")

DownloadData.R

This script can be used to download the phytoplankton and zooplankton time series data as well as the abiotic parameter.

PhenologyZP.R

In this script the Zooplankton metrics are calculated:

Imported raw abundance data from Sharkweb ¹
Biomass was calculated from individual bodymass based on location, life stage and sex according to the HELCOM COMBINE manual
Species were filtered to keep only the one with a good data coverage and excluding Nauplii
Density plot were performed based on the daily interpolation with a kernel smoothing
Weekly interpolation were made over the full time series with the Rpackage zoo and Week were retrieved using the package lubridate
Metric were calculated
- Peak timing using the center of gravity method
- Duration of the bloom using the cumulative biomass between the 25^th and 75^th percentile
- Magnitude of the bloom using the average of the biomass between the start and the end of the bloom
- A final dataset was saved combining all metric together “Output/Data/Weekly_Full_ZP.txt”

PhenologyPP.R

In this script the Phytoplankton metrics are calculated:

Imported raw biomass data from Sharkweb ²
Species were filtered to keep only the ones with a good data coverage and exclude heterotroph organisms
Density plot were performed based on the daily interpolation with a kernel smoothing
Spring bloom was the sum of the phytoplankton between Week 5 and 2
Summer bloom was the sum of the phytoplankton between Week 22 and 37
Weekly interpolation were made over the full time series with the Rpackage zoo and DOY were retrieved using the package lubridate
Metric were calculated
- Peak timing using the center of gravity method
- Duration of the bloom using the cumulative biomass between the 25^th and 75^th percentile
- Magnitude of the bloom using the average of the biomass between the start and the end of the bloom
A final dataset was saved combining all metric together “Output/Data/Weekly_Full_PP.txt”

Physical.R

In this script the Abiotic factor are calculated:

Imported raw data from Sharkweb ³
Weekly interpolation were made over the full time series with the Rpackage zoo and DOY were retrieved using the package lubridate
Monthly average are saved as “Output/Data/Abiotic.txt”

PlanktonSuccession.R

In this script the Fig. 1b and Fig. 3 are made:

Imported metrics and abiotic factor from
- “Output/Data/Weekly_Full_ZP.txt”
- “Output/Data/Weekly_Full_PP.txt”
- “Output/Data/Abiotic.txt”
Panel a shows the water temperature between 0 to 20m for each month at each station (average and range min-max)
Panel b shows the salinity between 0 to 20m for each month at each station (average and range min-max)
Panels c and d show total and taxa specific dynamics
- The peak timing (average between 2008 and 2022)
- The start and end of the bloom (average between 2008 and 2022)
A map is also plotted at the end of the script

Trends.R

In this script the Trends are analysed:

Imported data from
- “Output/Data/Weekly_Full_ZP.txt”
- “Output/Data/Weekly_Full_PP.txt”
- “Output/Data/Abiotic.txt”
Trends are analysed with the Kendall test from the Kendall package and effect size was assessed using the Sen’s slope from the mblm package for phytoplankton, zooplankton and seasonal salinity and temperature:
- Peak timing
- Peak initation
- Peak termination
- Peak magnitude
Coefficients of the trends for the phenological metrics are shown in Fig. 4 and Fig. S3
z-scores of the peak timing and magnitudes were calculated for all taxa and visualized across the year in Fig. 4 and Fig. S3
Coefficients of the trends and dynamics across the years for the abiotic parameters are shown in Fig. 2 and Fig. S2

GLM.R

In this script the Generalized Linear Models are made:

Imported data from
- “Output/Data/Weekly_Full_ZP.txt”
- “Output/Data/Weekly_Full_PP.txt”
- “Output/Data/Abiotic.txt”
Output of the GLM models are shown in Fig. 5 and Table S4

Nauplii.R

In this script Fig. S5 is produced showing seasonal abundance dynamics of Acartia, Pseudocalanus, Temora and Centropages nauplii stages are plotted

FreqPhytoProp.R

In this script the Fig. S1 and Fig. 1 are produced:

For all three stations, based on:
- Zooplankton sampling
- Phytoplankton sampling
- Temperature sampling
Density plot showing the total phytoplankton and zooplankton biomass dynamics and the contribution of each taxa

Validation_Trends.R

In this script Fig. S6 is produced:

We subsampled phytoplankton and zooplankton to have a similar temporal resolution than at BY15 and BY5
We weekly interpolated the data and calculated the phenology metrics for all taxa at BY31
We compared the subsampled data with the complete dataset

Data

First, the raw data were downloaded and are publicly available:

Manually by going on Sharkweb
Using the shaRk package and the script DownloadData

Note that data were downloaded using the shaRk Rpackage for this analysis. If the data are downloaded throughout the Sharkweb, make sure to update the headers accordingly.

The zooplankton bodymass were retrived from the HELCOM COMBINE manual, are shown in the Table S1 of the manuscript, and should have this structure to be compatible with the scripts:

#> 'data.frame':    132 obs. of  8 variables:
#>  $ Taxa   : chr  "Temora" "Temora" "Temora" "Temora" ...
#>  $ STAGE  : chr  "AD" "AD" "AD" "C1" ...
#>  $ SEXCO  : chr  "F" "M" "NS" "NS" ...
#>  $ Station: chr  "BY31" "BY31" "BY31" "BY31" ...
#>  $ Jan.Mar: num  60 55 57.5 9 18 60 55 57.5 9 18 ...
#>  $ Apr.Jun: num  65 40 52.5 6 15 65 55 60 6 15 ...
#>  $ Jul.Sep: num  60 45 52.5 6 15 40 45 42.5 6 15 ...
#>  $ Oct.Dec: num  50 50 50 4 14 50 45 47.5 6 15 ...

Where:

Taxa corresponds to the Taxa of interest
STAGE corresponds to the Stage of the organism (AD, C1, C4, JV or NS)
SEXCO corresponds to the sex code of the organism (F, M, NS)
Stationcorresponds to the station (BY31, BY15 or BY5)
Jan.Mar corresponds to the bodymass in Winter
Apr.Jun corresponds to the bodymass in Spring
Jul.Sep corresponds to the bodymass in Summer
Oct.Dec corresponds to the bodymass in Autumn

Processed data

In the folder Output/Data all the needed data can be found.

Weekly_Full_ZP.txt
Weekly_Full_PP.txt
Abiotic.txt

Weekly_Full_ZP.txt

This dataset is produced at the end of the script Phenology_ZP.R

It has this structure:

#> 'data.frame':    315 obs. of  11 variables:
#>  $ Taxa     : chr  "Acartia" "Acartia" "Acartia" "Acartia" ...
#>  $ Year     : int  2008 2008 2008 2009 2009 2009 2010 2010 2010 2011 ...
#>  $ Station  : chr  "BY15" "BY31" "BY5" "BY15" ...
#>  $ Group    : chr  "Copepoda" "Copepoda" "Copepoda" "Copepoda" ...
#>  $ Timing   : int  21 26 24 29 32 27 32 34 29 26 ...
#>  $ Maximum  : num  84.3 118.4 36.4 78.9 36.6 ...
#>  $ Magnitude: num  52 66.8 24.6 38 18.6 ...
#>  $ Start    : num  9 20 11.5 20 22.5 ...
#>  $ End      : num  28.3 33.5 35.5 34.8 40.7 ...
#>  $ Middle   : num  21 24.2 24.5 30 34 ...
#>  $ Duration : num  19.3 13.5 24 14.8 18.2 ...

Where:

Taxa: Mesozooplankton taxa: character
Year: Year between 2008 and 2022, integer
Station: the three stations BY31, BY15, and BY5, character
Group: Group of mesozooplankton taxa: Copepoda, Cladocera, and Rotatoria, character
Timing: Bloom peak timing (Week #), numeric
Maximum: Bloom magnitude maximum (\(\mu\)g L^-1), numeric
Magnitude: Bloom magnitude (\(\mu\)g L^-1), numeric
Start: Start of the growing season, correspond to the 25^th percentile of the cumulative biomass (Week #), numeric
End: End of the growing season, correspond to the 75^th percentile of the cumulative biomass (Week #), numeric
Middle: Middle of the growing season, correspond to the 50^th percentile of the cumulative biomass (Week #), numeric
Duration: Time span between the Start and the End of the growing season (Week #), numeric

Weekly_Full_PP.txt

This dataset is produced at the end of the script Phenology_PP.R

It has this structure:

#> 'data.frame':    404 obs. of  11 variables:
#>  $ Taxa     : chr  "Cyanobacteria" "Cyanobacteria" "Cyanobacteria" "Dinoflagellates" ...
#>  $ Year     : int  2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
#>  $ Station  : chr  "BY15" "BY31" "BY5" "BY15" ...
#>  $ Group    : chr  "Phytoplankton" "Phytoplankton" "Phytoplankton" "Phytoplankton" ...
#>  $ Timing   : int  26 31 30 17 17 19 25 26 24 21 ...
#>  $ Maximum  : num  60.3 14.3 117 67.7 43.4 ...
#>  $ Magnitude: num  55.66 8.44 91.21 59.57 30.2 ...
#>  $ Start    : num  23.7 24.8 27.6 14.5 14.8 ...
#>  $ End      : num  29 34.2 30.9 19.1 18 ...
#>  $ Middle   : num  26.4 29.8 29.5 16.7 16.6 ...
#>  $ Duration : num  5.33 9.45 3.25 4.6 3.23 ...

Where:

Taxa: Phytoplankton taxa, See Table S3 for more information character
Year: Year between 2008 and 2021, integer
Station: the three stations BY31, BY15, and BY5, character
Group: Group of phytoplankton taxa: only Phytoplankton, character
Timing: Bloom peak timing (Week #), numeric
Maximum: Bloom magnitude maximum (\(\mu\)gC L^-1), numeric
Magnitude: Bloom magnitude (\(\mu\)gC L^-1), numeric
Start: Start of the growing season, correspond to the 25^th percentile of the cumulative biomass (Week #), numeric
End: End of the growing season, correspond to the 75^th percentile of the cumulative biomass (Week #), numeric
Middle: Middle of the growing season, correspond to the 50^th percentile of the cumulative biomass (Week #), numeric
Duration: Time span between the Start and the End of the growing season (Week #), numeric

Abiotic.txt

This dataset is produced with the script Physical.R

It has this structure:

#> 'data.frame':    4600 obs. of  6 variables:
#>  $ Year     : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
#>  $ month    : int  1 1 1 1 1 1 1 1 1 1 ...
#>  $ Parameter: chr  "Nutrient" "Nutrient" "Nutrient" "Oxy_0.20" ...
#>  $ Station  : chr  "BY15" "BY31" "BY5" "BY15" ...
#>  $ Value    : num  2.64 3.81 3.04 7.92 7.85 ...
#>  $ n        : int  1 1 1 1 1 1 1 1 1 1 ...

Where:

Year: Year of the sampling event, integer
month: Month of the sampling event. 1 corresponds to January, 2 to February, and so on, integer
Parameter: Parameter sampled. character
- Oxy_0.20: Integrated oxygen concentration between 0 and 20m depth, mL L^-1
- SST: Sea surface temperature, °C
- Temp_0.20: Integrated temperature between 0 and 20m depth, °C
- Salinity: Integrated salinity between 0 and 20m depth
- Salinity_SSS: Sea surface salinity
- Nutrient: Nitrate and nitrite concentration between 0 and 20m depth, \(\mu\)mol L^-1
- strat_index: Stratification index calculated as (SST - Temperature at 60m) / average temperature between 0 and 60m
Station: Station sampled, character
Value: value of the parameter, numeric
n: number of observation, integer

Note that data were downloaded using the shaRk Rpackage for this analysis. If the data are downloaded throughout the Sharkweb, make sure to update the headers accordingly.↩︎
Note that data were downloaded using the shaRk Rpackage for this analysis. If the data are downloaded throughout the Sharkweb, make sure to update the headers accordingly.↩︎
Note that data were downloaded using the shaRk Rpackage for this analysis. If the data are downloaded throughout the Sharkweb, make sure to update the headers accordingly.↩︎